clc;
clear all;

%%%%% LOAD DATA %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    % load fips codes
    fips = csvread('../../data/tradecosts/Fips.csv')';

    % load fips codes
    avgship = csvread('../../data/tradeflows/avg_value_per_shipment.csv')';

    % set number of counties
    Counties = length(fips);

    % load inflation data
    P = csvread('../../raw/tradecosts/tradecost_inflation.csv');

    % load fuel cost data
    F = csvread('../../raw/tradecosts/tradecost_fuelcosts.csv');
    C = F(:,1).*(1./F(:,2)).*P;

    % load wage cost data 
    W0 = csvread('../../raw/tradecosts/tradecost_wagecosts.csv');
    W = W0.*P;

    % load income data
    Y0_mfg = csvread(['../../raw/tradecosts/Y_mfg.csv']);
    Y0_rem = csvread(['../../raw/tradecosts/Y_rem.csv']);
    
    % load labor data
    L0_mfg = csvread(['../../raw/tradecosts/L_mfg.csv']);
    L0_rem = csvread(['../../raw/tradecosts/L_rem.csv']);
 
    % load trade cost data
    Dist1980 = csvread('../../data/tradecosts/Distance1980.csv');
    Time1980 = csvread('../../data/tradecosts/Time1980.csv');
    Dist1990 = csvread('../../data/tradecosts/Distance1990.csv');
    Time1990 = csvread('../../data/tradecosts/Time1990.csv');
    Dist2000 = csvread('../../data/tradecosts/Distance2000.csv');
    Time2000 = csvread('../../data/tradecosts/Time2000.csv');
    
    % choose theta 
    theta = 4;
    
%%%%% LOOP OVER ALL YEARs %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    for t = 1980:2000

        yearFact = num2str(t);
        disp(['The year is ',num2str(t)]);
        if t >=1969

            if t < 1975

                % load trade cost data, by year and cost parameters
                Time = zeros(max(size(fips)),max(size(fips)));
                Dist = zeros(max(size(fips)),max(size(fips)));
                DistFact = Dist1970;
                for i = 1:max(size(fips))
                Dist(1:max(size(fips)),i) = DistFact(1+(i-1)*max(size(fips)):i*max(size(fips)));
                end
                DistFact = Dist;
                DistCost = DistFact.*C(t-1979,1);
                TimeFact = Time1970;
                for i = 1:max(size(fips))
                Time(1:max(size(fips)),i) = TimeFact(1+(i-1)*max(size(fips)):i*max(size(fips)));
                end
                TimeFact = Time;
                TimeCost = TimeFact.*W(t-1979,1);
                tauFact  = ((DistCost + TimeCost));
                
            end

            if t >= 1975 

                if t <= 1985

                    % load trade cost data, by year and cost parameters
                    Time = zeros(max(size(fips)),max(size(fips)));
                    Dist = zeros(max(size(fips)),max(size(fips)));
                    DistFact = Dist1980;
                    for i = 1:max(size(fips))
                    Dist(1:max(size(fips)),i) = DistFact(1+(i-1)*max(size(fips)):i*max(size(fips)));
                    end
                    DistFact = Dist;
                    DistCost = DistFact.*C(t-1979,1);
                    TimeFact = Time1980;
                    for i = 1:max(size(fips))
                    Time(1:max(size(fips)),i) = TimeFact(1+(i-1)*max(size(fips)):i*max(size(fips)));
                    end
                    TimeFact = Time;
                    TimeCost = TimeFact.*W(t-1979,1);
                    tauFact  = ((DistCost + TimeCost));
                 
                end

                if t >= 1986

                    if t <= 1995

                        % load trade cost data, by year and cost parameters
                        Time = zeros(max(size(fips)),max(size(fips)));
                        Dist = zeros(max(size(fips)),max(size(fips)));
                        DistFact = Dist1990;
                        for i = 1:max(size(fips))
                        Dist(1:max(size(fips)),i) = DistFact(1+(i-1)*max(size(fips)):i*max(size(fips)));
                        end
                        DistFact = Dist;
                        DistCost = DistFact.*C(t-1979,1);
                        TimeFact = Time1990;
                        for i = 1:max(size(fips))
                        Time(1:max(size(fips)),i) = TimeFact(1+(i-1)*max(size(fips)):i*max(size(fips)));
                        end
                        TimeFact = Time;
                        TimeCost = TimeFact.*W(t-1979,1);
                        tauFact  = ((DistCost + TimeCost));
                       
                    end

                    if t >= 1996

                        if t <= 2005 

                            % load trade cost data, by year and cost parameters
                            Time = zeros(max(size(fips)),max(size(fips)));
                            Dist = zeros(max(size(fips)),max(size(fips)));
                            DistFact = Dist2000;
                            for i = 1:max(size(fips))
                            Dist(1:max(size(fips)),i) = DistFact(1+(i-1)*max(size(fips)):i*max(size(fips)));
                            end
                            DistFact = Dist;
                            DistCost = DistFact.*C(t-1979,1);
                            TimeFact = Time2000;
                            for i = 1:max(size(fips))
                            Time(1:max(size(fips)),i) = TimeFact(1+(i-1)*max(size(fips)):i*max(size(fips)));
                            end
                            TimeFact = Time;
                            TimeCost = TimeFact.*W(t-1979,1);
                            tauFact  = ((DistCost + TimeCost));
                          
                        end

                        if t >= 2006 

                            % load trade cost data, by year and cost parameters
                            Time = zeros(max(size(fips)),max(size(fips)));
                            Dist = zeros(max(size(fips)),max(size(fips)));
                            DistFact = Dist2010;
                            for i = 1:max(size(fips))
                            Dist(1:max(size(fips)),i) = DistFact(1+(i-1)*max(size(fips)):i*max(size(fips)));
                            end
                            DistFact = Dist;
                            DistCost = DistFact.*C(t-1979,1);
                            TimeFact = Time2010;
                            for i = 1:max(size(fips))
                                Time(1:max(size(fips)),i) = TimeFact(1+(i-1)*max(size(fips)):i*max(size(fips)));
                            end
                            TimeFact = Time;
                            TimeCost = TimeFact.*W(t-1979,1);
                            tauFact  = ((DistCost + TimeCost));
                           
                        end

                    end 

                end

            end

        end
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        % Tau 
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

        % save tau for this year
        tauFact = 1 + (tauFact./avgship(1,1));
        dlmwrite(['../../data/tradecosts/Tau' yearFact '.csv'],[fips', tauFact],'delimiter',',','precision',17);

        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        % Manufacturing 
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        % load income for this year
        Y_mfg = Y0_mfg(:,t-1979)';

        % choose labor for this year
        L_mfg = L0_mfg(:,t-1979)';

        % solve for market access for this year 
        maFact_mfg = SolveForMA(L_mfg, tauFact, theta, Counties);

        % save market access so far
        MA_mfg(:, t - 1979) = maFact_mfg;
        dlmwrite(['../../data/tradecosts/MA_mfg.csv'],[fips', MA_mfg],'delimiter',',','precision',17);
        
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        % Manufacturing 
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        % load income for this year
        Y_rem = Y0_rem(:,t-1979)';

        % choose labor for this year
        L_rem = L0_rem(:,t-1979)';

        % solve for market access for this year 
        maFact_rem = SolveForMA(L_rem, tauFact, theta, Counties);

        % save market access so far
        MA_rem(:, t - 1979) = maFact_rem;
        dlmwrite(['../../data/tradecosts/MA_rem.csv'],[fips', MA_rem],'delimiter',',','precision',17);
       
    end